Quantum Dimer Model on the triangular lattice: Semiclassical and variational 
approaches to vison dispersion and condensation 



Gregoire Misguich* 
Institut de Physique Theorique, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France 

Frederic Mila 

Institute of Theoretical Physics, Ecole Polytechnique Federate de Lausanne, CH-1015 Lausanne, Switzerland 

(Dated: January 23rd, 2008) 

After reviewing the concept of vison excitations in Z2 dimer liquids, we study the liquid-crystal 
transition of the Quantum Dimer Model on the triangular lattice by means of a semiclassical spin- 
wave approximation to the dispersion of visons in the context of a "soft-dimer" version of the model. 
This approach captures some important qualitative features of the transition: continuous nature of 
the transition, linear dispersion at the critical point, and vl2 x \/T2 symmetry-breaking pattern. In 
a second part, we present a variational calculation of the vison dispersion relation at the RK point 
which reproduces the qualitative shape of the dispersion relation and the order of magnitude of the 
gap. This approach provides a simple but reliable approximation of the vison wave functions at the 
RK point. 
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I. INTRODUCTION 

Since they have been shown to possess Resonating Va- 
lence Bond (RVB) phases on the triangular, 1 kagome 
and other (non-bipartite) lattices, ' Quantum Dimer 
Models (QDM) have been one of the main paradigms in 
the field of quantum spin liquids. These models, where 
the Hilbert space is spanned by hard-core dimer coverings 
of the lattice, are expected to capture the phenomenology 
of quantum antiferromagnets where the wave function is 
dominated by short-range valence bond configurations. 
On the triangular lattice, the simplest QDM is defined 
by the Hamiltonian: 

"= -< £(Lv)Crl +H - c -) 
+"£(l//)(//H.r)Crl)< i > 
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where the sum runs over all plaquettes (rhombi) includ- 
ing the three possible orientations. The kinetic term, con- 
trolled by t, flips the two dimers on every flippable pla- 
quette, i.e., on every plaquette with two parallel dimers, 
while the potential term controlled by the interaction V 
describes a repulsion (V > 0) or an attraction (V < 0) 
between nearest-neighbor dimers. 

The RVB phase is now relatively well understood. The 
first result goes back to Rokhsar and Kivelson who 
showed that, for V/t = 1 (the Rokhsar-Kivelson or RK 
point), the ground state is the sum of all configurations 
with equal amplitudes: 

i^> = ^Ei c >- ( 2 ) 

Since then, dimer-dimer correlations have been shown to 
be short-ranged in a range of parameters below V/t — 



l, 1 ' 5 and the excitation spectrum to be gapped. 1 ' 1 '''' 8 Be- 
sides, it has topologically degenerate ground states on 
non-simply connected clusters. '°'' This degeneracy is 
not related to a standard symmetry breaking. Indeed, 
there is no local order parameter, 1 '" but only Z 2 topo- 
logical order. 

In QDM, the nature of the phase transition from a liq- 
uid to a solid is a long standing problem. It goes back to 
Jalabert and Sachdev 11 (see also 1 ") who studied a three- 
dimensional frustrated Ising model related to the square 
lattice QDM. 

In a more general context, Senthil and Fisher 1,5 ' 14 
showed that models of Mott insulators can be cast in the 
form of a Z 2 gauge theory. In this language, the tran- 
sition from a fractionalized insulator to a conventionally 
ordered insulator appears to be a condensation of Z 2 - 
vortices (dubbed visons). Using a duality relation, they 
showed that such transitions correspond to an ordering 
transition in a frustrated Ising model in transverse field. 
As we will see, this applies to the present QDM. 

Building on a mapping between QDM's at V = and 
Ising models in a transverse field, Moessner et al. le ' 
have developed a Landau-Ginzburg approach and sug- 
gested that the transition could be continuous and in the 
three-dimensional 0(4) universality class. 

Numerical evidence in favor of this scenario has been 
obtained with Green's function Quantum Monte Carlo 
by Ralko et al., who have shown that, at the transition 
point, the static form factor of the crystal decreases to 
zero on the crystal side of the transition, while the dimer 
gap also decreases to zero on the liquid side. 18 More re- 
cently, the vison spectrum has also been numerically de- 
termined using Green's function Quantum Monte Carlo, 
with the conclusion that a soft mode indeed develops at 
the transition. 

In spite of these results, a simple picture for the wave 
functions of the visons and the evolution of their spec- 



trum is still missing. One difficulty is that, like any 
vortex, these excitations cannot be created by operators 
which are local (in the dimer variables). To attack this 
problem, we follow two strategies. First of all, we look at 
the problem in the context of a Z2 gauge theory on the 
triangular lattice. As usual, this theory can be mapped 
onto a dual Ising model. The duality transforms the 
nonlocal excitations of the gauge theory into local exci- 
tations, which we study using a semiclassical approach. 
This simple 1/5 expansion already captures most as- 
pects of the confinemcnt-dcconfincment transition (soft 
mode and condensation) of the Z2 gauge theory. Sec- 
ondly, building on the explicit form of the excitations 
of the Z2 gauge theory, we construct single vison wave 
functions for the QDM and study the properties of the 
system at the RK point in the variational Hilbert space 
spanned by these states. These "variational visons" - liv- 
ing on the triangular plaquettes and experiencing a flux 
7r emanating from each site of the lattice - turn out to 
be linearly independent. The associated dispersion is in 
good qualitative agreement with that obtained in Monte 
Carlo simulations and the value of the gap at the RK 
point (A var = 0.119), obtained with a single variational 
parameter, has the correct order of magnitude (Monte 
Carlo simulations'' give 0.089). Improving quantitatively 
further these variational results would require more ad- 
justable parameters to account in more details for the lo- 
cal (dimer-dimer, etc) correlations in the vicinity of the 
core of the vortex, a task which has not been carried out 
here. 



II. VISONS IN Z 2 GAUGE THEORY 



The connection between QDMs and Z2 gauge theories 
has already been discussed from different perspectives 
(see in particular Refs. 2,3,12,20), and is rooted in the 
existence, in both families of models, of Ising-like degrees 
of freedom subjected to local constraints (hard-core con- 
straints for the dimers, Gauss law for the gauge theory). 

In this section, we review two known mappings: i) 
from a Z 2 lattice gauge theory on the hexagonal lat- 
tice to the triangular lattice QDM (valid in a particular 
limit, Sec. II A) and ii) from the Z2 theory to its dual 
frustrated Ising model (Sec.IIB). The latter Ising model 
is then studied using a semiclassical approximation in 
Sees. II B 2 and II B 3. Since the Z 2 gauge theory - QDM 
mapping is formally justified only in the confined phase 
of the gauge theory, the relevance of the results obtained 
in the RVB (deconfined) phase of the QDM is not guar- 
anteed a priori. Some reasons why this is expected to 
be the case are discussed in the next section, when we 
build on the results obtained for the excitations of the 
Z 2 gauge theory in its deconfined phase to construct el- 
ementary excitations in the RVB phase of the QDM. 
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FIG. 1: (Color online) Some useful definitions of sets of bonds: 
1) is the set of the three bonds forming the edges of the 
plaquette i. 2) l[a] is the set of the six (fat) bonds emanating 
from the triangular site a. 3) I < i denotes the bonds forming 
a "zigzag" string extending (to the left) from the triangular 
plaquette i to an edge of the lattice. 



A. Z2 gauge theory 

To write down a Z2 gauge theory analogous to a QDM, 
the starting point is to define Pauli matrices f; on the 
bonds I of the triangular lattice such that r'f = —1 if 
bond I is occupied by a dimer and +1 if it is empty, 
while rf changes the state of bond I. Since each kinetic 
term of the QDM flips the dimers around a rhombus, it 
corresponds to a product of four rf around this rhombus. 
This would, however, not keep the structure of the simple 
Z2 gauge theories, in which the kinetic term acts on an 
elementary plaquette of the lattice, a crucial ingredient 
to get a simple Ising model by duality. An alternative 
is to consider the more standard Z 2 gauge theory which, 
in its Hamiltonian formulation, is defined by: 



H = Hj + H T 



-'£*f- r £ir 



i(z)> 



(3) 



where i runs over the sites of the dual honeycomb lat- 
tice, and l(i) are the three bonds forming the triangular 
plaquette around site i (see Fig. 1). The hallmark of this 
model is to have local conserved quantities. Indeed, 



J [a] J 



0. 



(4) 



i[o] 



where a is a site of the triangular lattice, and the prod- 
uct over / [a] runs over the six links emanating from a (see 
Fig. 1). This allows one to define different sectors accord- 
ing to whether Yli[ a ] T [ a \ ^ s ec l ua l to +1 or — 1. Since in 
the QDM the number of dimers emanating from a given 
site is exactly equal to 1, it is clearly better to consider 
the sector where 



n 



= -1 



(5) 



for all a since this forces the number of dimers emanating 
from a site to be odd (defining an odd Ising gauge theory 
in the terminology of Ref. ). Then, the true constraint 
is recovered in the limit T/J — > if J > since, in the 
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ground state, the number of dimers is then minimal. A 
bona fide QDM is then recovered if Hj> is treated with 
degenerate perturbation theory. Since H r changes the 
number of dimers, its effect vanishes to first order. To 
second order however, one recovers exactly the QDM of 
Eq. (1) with V = and t = T 2 /J. So, if J/T > 1, the Z 2 
gauge theory maps onto the QDM at V/t = 0. At finite 
J/T the model can be viewed as a "soft-dimer" model, 
where 1, 3 or 5 dimers may touch a given site. 

As we shall see below, the limit J/T 3> 1 lies deep 
inside the confined phase of the gauge theory, and from 
previous work on the QDM, it is known that at V/t = 0, 
the model is in a Valence Bond Crystal phase. It would, 
of course, be very interesting to connect the two models 
away from this limit, when the gauge theory is in its 
deconfmed phase and the QDM is in the dimer liquid 
phase. An interesting step in this direction is provided 
by higher order perturbation theory. Indeed, to fourth 
order in T/J, a repulsion between dimers V — T 4 /2J 3 
is generated. However, other terms of the same order, 
involving dimer shifts along loops of length 6, are also 
generated (see Appendix A). So, a rigorous mapping does 
not extend beyond the V/t — point. 




FIG. 2: (Color online) Dashed lines: Horizontal paths used 
in Eq. 6 connecting the boundary of the system to triangles 
i and j. Dotted segment: Path Q^. In this example 
crosses one dimer of the reference configuration Co (ellipses) 
and, therefore, ey = — 1 (Eq. 48) 

leading to the same Ising model up to a gauge transfor- 
mation. The present choice leads to the smallest unit cell 
(4 honeycomb sites). 

In the following, we study the phase diagram and the 
excitations of this model within a semiclassical (large S) 
approximation. 



B. Dual Ising model 

As usual, this %i gauge theory is best analyzed by 
mapping it onto a dual Ising model in a transverse field. 

1. The model 

This can be achieved by introducing spin-^ operators 
cji on the dual honeycomb lattice: 

^ = Ihw> < = ( 6 ) 

l(i) Ki 

where / < i represents all bonds cutting a straight path 
(say, horizontal, see Fig. 2 or 1) starting at i (we im- 
plicitly assume a finite lattice with open boundary con- 
ditions). Combined with the constraint, this definition 
implies that 

o-i„] = M l3 7f (7) 

for two neighboring sites i,j separated by the bond 
Mij = ±1 is such that each hexagon has exactly one 
Mij = — 1 bond (see Fig. 3). In terms of these spin 
operators, the Hamiltonian is the fully frustrated Ising 
model (FFIM) discussed by Moessner and Sondhi: 1 '' 26 

H = Hj + H r = - J Mijo-fa* -T^af. (8) 

{id) 

Note that choosing other paths to define of leads to other 
signs for Mj,-, which, however, are always such that an 
odd number of minus signs appear around each hexagon, 



2. Classical phase diagram 

First, we determine the classical ground state of the 
model as a function of T/J. To this end, we replace 
the spin-1 operators by classical 3-component vectors of 
unit length. Since the y component does not appear in 
the Hamiltonian, it is clear that in the ground state the 
magnetization must lie in the x-z plane. To find the low- 
est energy solution, we use the following parametrization: 

of - P(i), (9) 
of - VI-pW 2 , (10) 

with \p(i) \ < 1. The corresponding energy is 

E = -\ E M ijp(Mj) - r£ Vi^pW- (ii) 

ij i 

Initially defined for nearest neighbors (Eq. 7) , the coef- 
ficients My have been upgraded to a matrix M by setting 
all other elements to 0. This matrix describes the mo- 
tion of a particle on a honeycomb lattice with 4 sites per 
unit cell and a flux 7r per hexagonal plaquette. It re- 
duces to a 4x4 matrix after Fourier transformation. The 
eigenvalues associated with the momentum k = (k x ,k y ) 
are 17 



± y 3 ± y/2[3 + cos^fcj;) — cos(k x + k y ) + cos(k x — k y )\. 

Likewise, we denote by p the column vector of compo- 
nents {pi}. 
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FIG. 3: (Color online) The fat (green) bonds of the hexagonal 
lattice (crossed by a dimer of the reference configuration on 
the triangular lattice) have My = — 1, the other bonds have 
Mij = 1. The four sites of the unit cell (small blue triangles) 
are labelled from 1 to 4 and the Bravais vector are x and y. 



T/J > v6 For r >• J, we expect the p(i)'s to be 
We can, therefore, expand Eq. 11 to 



small (or zero), 
quadratic order 



E 



J 



p l Mp- 



1 



(12) 



The largest eigenvalue of M being we find that, as 
long as T/J > y/6, the energy is minimized by p = and 
all spins point in the x direction. 

< T/J < \f% At T/ J = y/6, all the real eigenvectors 

of M for the eigenvalue \/6 (satisfying \pi\ < 1) minimize 
Eq. 12. Let us choose four complex vectors which form 
a basis of the subspace associated with the eigenvalue 
V6: 17 



vi(x } y) = v* 3 



V2{x,y) = v* 4 = 



F = 



e 5i7r/12 ip 
1 

e -ijr/12 

e i7r/12 
5i7r/6 

1/F 

e -U-K/12 lp 

2sin(57r/12), 



exp 



exp 



in it 
~6 X + l 2 y 



5ir 7r 

i— x + i—y 



(13) 



where the 4 entries of the vector refer to the 4 sites in 
the unit cell (numbered as in Fig. 3) and x,y are the 
(Bravais) coordinates of the unit cell. These eigenvectors 
correspond to the four points labelled B in the rectangu- 
lar Brillouin zone of Fig. 6. 

The most general real eigenvector p can be 
parametrized by three angles ai, ct2, and (3 and a nor- 
malization factor A: 

p{i) = A [cos(f3)n(v 1 (i)e lai ) + sm(l3)n(v 2 (i)e la2 )] . 

(14) 

To find the ground state, the energy has to be minimized 
as a function of these 3 parameters. The analysis is made 
easier when one realizes that the second and fourth mo- 
ments of the spin deviations are independent of the three 



angles. Indeed, 

i 



with r 2 



r 2 A 2 , 
2r 4 A 4 , 



1 f 1 
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(15) 
(16) 
(17) 



Replacing p by Eq. 14 in the expression for the energy 
(Eq. 11), we get up to order A 4 : 

E/N = -T+^(r - jV6jr 2 X 2 + ^Tr i X 4 + 0(X 6 ) (18) 



Minimizing with respect to A 2 gives: 



A 2 = 



/ 2 r ' 



(19) 



For r > the energy is minimized for A = 0. For 

r < we have to expand Eq. 11 to the 6 th order to 

find the angles ai, «2 and (3 which minimize the classical 
energy. Indeed, Eq. 18 shows that, up to order A 4 , the 
energy is independent of the three angles. At this order, 
the energy is constant and minimum on a 3-dimensional 
sphere, a consequence of the 0(4) symmetry discovered 
by Moessner and Sondhi. One, therefore, has to go to 
the next order to find the actual minima. At 6 th order 
in A, the energy is minimized when 



(20) 



is minimum. A numerical investigation shows that the 
solutions (for T/J close to but below V6) are 48-fold 
degenerate and can be deduced from each other by sym- 
metry operations (12 translations and 4 point-group op- 
erations). This confirms the result obtained previously 
on symmetry grounds. 11 '' 17 

Motivated by the FFIM/dimer model correspondence, 
we are interested in the average "dimer density" : 



(21) 



for all pairs of nearest-neighbor honeycomb sites. 

In the classical limit, this may be approximated by 



1 



(l-M ijP (i)p(j)), 



(22) 



where p (Eq. 14) is the classical ground state. Close to 
the transition at T c = J\/6, p scales as A ~ ^/T c — T 
(Eq. 19) and the "dimer density" , thus, shows small de- 
viations about tj. To visualize the "dimerization" pat- 
tern in the vicinity of r c , the appropriate quantity is, 
therefore, a relative rescaled "dimer" density defined by 



D 4 



ATI 



(dij — i), and plotted in Fig. 4 for one of the 
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FIG. 4: (Color online) Plot of the rescaled relative "dimer" 
density Dij in one of the 48 classical ground states of the 
FFIM in transverse field, for a transverse field V just below 
F c . Thin red bonds represent Dij — (corresponding to the 
highest "dimer density" dij = 0.5). Blue bonds are for Dij < 
0, that is, a lower dimer density, with a width proportional 
to \Dij\. Dij takes only 4 different values: 0, -0.259, -0.518, 
and 



-0.776. 



ground states. The obtained pattern is highly reminis- 
cent of the VT2 x \/l2 VBC observed in the triangular 
lattice QDM. In particular, the 48 spin configurations 
give only 12 different "dimer" patterns, with a smaller 
unit cell containing 12 sites of the triangular lattice. It 
is also interesting to notice that p(i) vanishes in the tri- 
angles located inside the large diamonds (marked with 
a dot in Fig. 4). From Eq. 10, the corresponding sites 
have a magnetization pointing exactly in the x direction, 
which means a total absence of "vison" and local correla- 
tions identical to that of the J = paramagnetic ground 
state. The "dimer" density dij is uniform inside these 
diamonds, in qualitative agreement with the intriguing 
observation made in Ref. 18. 



3. Semi- classical analysis for T/J > \f& 

If J = 0, the ground state is the fully polarized state 
af = 1. The first (degenerate) excited state is obtained 
by flipping one spin at some arbitrary triangle. To first 
order in J/T, the spectrum is obtained by diagonalizing 
the perturbation Hj in this subspace. This amounts to 
diagonalizing M, and leads to the eigenvalues e k — 2T + 
Jm l k , a result already obtained previously. 

To go further, beyond the limit J/T <C 1, we perform a 
1/S semiclassical expansion. We generalize the spin-1/2 
Hamiltonian to an arbitrary value S of the spin: 



^ (W) 



(23) 



(which reduces to Eq. 8 when S — 1/2). Spin devia- 



tions away from the x directions are represented using 
Holstein-Primakoff bosons. To leading order in 1/5: 



1 



s 



btbi 



(24) 
(25) 



One truncates the Hamiltonian to quadratic order in bi, 
and diagonalizes it through a Bogoliubov transformation: 



(26) 



For the a - operators to be bosonic creation operators, 

the matrices U and V must satisfy Wll — V^V = 1. 
It is convenient to introduce a unitary matrix f2 which 
transforms M into a diagonal matrix M: 



M = VLMn\ 
M kp = m k 5k P - 



(27) 
(28) 



In this new basis, we can look for diagonal solutions for 
U and V: 



u = nurt , v = nvvt\ 

Uk P = UkSkp , Vk p — VkSkp- 



(29) 
(30) 



Since M is real, we can choose S O(N). After some 
algebra, one finds that the u k and v k which diagonalize 
H are given by 



tanh(26> fe ) 



cosh(6' fe ), 
sinh(6» fc ), 



(31) 
(32) 

(33) 



r^. + 2r/j 

The energies of the Bogoliubov excitations are given by 



ffc 




m k 



(34) 



As expected, the spectrum becomes gapless at T/J = v6 
and we recover the localized spin flip energy e k ~ T/S 
when r — > oo. The Bogoliubov spectrum is shown Fig. 5 
for a few values of T at or above r c . One clearly sees 
that the dispersion is linear around the gapless points at 
the transition. 

This behavior is remarkably similar to that found nu- 
merically by Ralko et al. V) for the QDM by Green's func- 
tion Quantum Monte Carlo. In the next section, we build 
on this resemblance to develop a variational approach to 
the vison spectrum of the QDM. Note that, as stated 
above, the Z 2 gauge theory and the QDM can only be 
rigorously mapped onto each other deep into the con- 
fined (resp. VBC) phase, and the resemblance between 
their spectra in the deconfined (resp. RVB) phases might 
seem fortuitous. A somewhat deeper connection will be 
described in the next section. 
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FIG. 5: (Color online) Spin- wave dispersion relation Set/ J 
(Eq. 34) for the path A^B^C^A in the Brillouin zone (see 
Fig. 6) and different values of T/J:3 (top), 2.75, 2.5, and 
2d 2.448 (bottom). At the critical point T/J = \/6, the 
spectrum is linear around k — (it/6, 7r/2) and k = (57r/6, 7r/2). 
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FIG. 6: (Color online) Dashed rectangle: Brillouin zone 
(— 7r < k x < 7r and — n < k y < n) of the hexagonal lattice 
shown in Fig. 3. A, B, and C are the high-symmetry points 
used in Fig. 5. The large hexagon is the Brillouin zone of the 
underlying triangular lattice. 



III. VISONS IN THE QUANTUM DIMER 
MODEL 

In the Ising model, elementary excitations for large 
enough T/J are spin flips, induced by erf , which delocal- 
ize and get dressed under the effect of H j . In the equiv- 
alent dual gauge theory, this excitation is produced by 
the nonlocal string operator of = n;<i T f ■ By analogy, 
it is natural to define a "point" vison creation operator 
for the QDM by 5 ' 22 

Vi = (-lf^ , ft,, (35) 

Ki 

where the dimer operator h\ is defined by hi = 1 if bond 
I is occupied, and otherwise. Again, we consider here 
a finite lattice with open boundary conditions. Let us 
first see to which extent this operator is the analog of 
the vortex creation operator in Ising gauge theories. 



A. Z2 gauge structure of QDM 

In gauge theories, the Wilson loop operator defined by 

Wan = I] T i ( 36 ) 
lean 

plays a central role since it allows one to distinguish (in 
the absence of matter field, as here) the deconfined and 
confined phases depending on whether its ground state 
expectation value (or flux) tends to zero exponentially 
with the perimeter of the domain £1 (deconfined) or with 
the area of the domain (confined). For a Z 2 gauge theory 
Wqq = 1 and the flux going through Q measured by Won 
can only take two values ±1. In that respect, an inter- 
esting property of the operator of = Jlz<i T i YS that it 
changes the flux between -1 and +1 if the site i is inside 
the domain. Then, of oj commutes with the Wilson loop 
and does not change the flux unless i and j sit on op- 
posite sides of the boundary. Since the deconfined phase 
can be accessed from the J<T limit, the ground state 
can be thought of as an Ising paramagnet (cr x = 1 and 
Wgn = 1 everywhere) "dressed" perturbatively by suc- 
cessive applications of Jo~fo~j, where i and j are nearest 
neighbors. Since Jof cr| only changes the flux for pairs ij 
across the boundary, the expectation value of the Wilson 
loop operator behaves according to a perimeter law. A 
similar pcrturbative argument can be used to derive the 
area law in the confined phase. 

Most of these standard Z 2 gauge theory results apply 
to the QDM, with one difficulty however. The Wilson 
loop operator cannot be defined in the same way since 
flipping the dimer occupation along a loop will often 
lead to an unphysical state that violates the condition 
of having exactly one dimer emanating from each site. 
Different ways to overcome this problem can be envis- 
aged. One possibility is to accept that the Wilson loop 
operator gives zero when applied to states that are non- 
flippable along the loop. The property Wg n = 1 is lost, 
but Wgn has eigenvalues 0, +1, or —1, and eigenstates 
with eigenvalues +1 or —1 are still interchanged when a 
vison operator (Eq. 35) is applied inside the domain. So 
confinement or deconfinement of visons is still expected 
to lead to area or perimeter laws. 

Alternatively, starting from the observation that, when 
it does not give zero, the Wilson loop operator just 
shifts the dimers along the contour, one can try to de- 
fine a (complicated) flux operator in the dimer space 
which shifts the dimer along a "fattened" contour which 
"adapts" locally to the dimer configuration it acts upon. 
Although tricky to define in practice,"' this operator 
would have the advantage that its square is still equal 
to 1 , preserving the manifest Z 2 structure of the theory. 

Another way to underline the deep connection between 
the two models is to consider the Ising model as a "soft- 
dimer" model and introduce the projection operator onto 
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the "hard-core" dimer Hilbert space, defined by: 

p = TT (n p -3)(n p -5) 
11 (l-3)(l-5) 



(37) 



where the operator n p counts the number of frustrated 
bonds (dimers) around the plaquette p: 



1 6 

-^(1-M iii+1 (7?< +1 ) 



(38) 



By construction, P = 1 on all the Ising configurations 
which correspond to a valid "hard-core" dimer covering 
(of the triangular lattice) and P = otherwise. Now, P 
commutes with af since all terms in P contain an even 
number of a z operators. This means that P conserves 
the total flux, so that a spin state with -1 flux becomes a 
dimer wave function with an odd number of visons after 
projection. 

B. Point vison 

Equation 35 defines the simplest operator which 
changes the flux, that is, which creates a vison. So we 
may consider 

= V{i)\RK) = 4ftY,(- l ) m \ c ) ( 39 ) 



Af ■ 



1 



Ar?|{of = +1}) 



(40) 



as a first variational approximation to the true lowest 
eigenstate of H. in the —1 flux sector (Af is the total 
number of dimer coverings). The ground state \RK) is 
a zero-energy eigenstate of TL, implying that the expec- 
tation value Ek < of the kinetic energy term exactly 
compensates the expectation value Ep > of the poten- 
tial energy term. Indeed, 



E F 



3N Po 



-E 



K, 



(41) 



where 3 N is the total number of rhombi (iV the number of 
sites) and po the probability to have 2 parallel dimers on 
a given rhombus in the classical dimer problem (with uni- 
form measure over all dimer configurations). Using the 
Pfaffians method, ' one finds in the thermodynamic 
limit 

Po ^ 0.0933310104... (42) 

In |i) the expectation value of the potential energy 
term is the same as in \RK). In fact, all the observables 
which are diagonal in the dimer basis commute with a z 
and thus have the same expectation value in the ground 
state \RK) and in the trial wave function The in- 
crease of the energy is only kinetic. Because the kinetic 
terms corresponding to the three diamonds around i an- 
ticommute with af , their expectation values has changed 
sign. Thus we find 



(i\H\i) - (RK\H\RK) = 6p ~ 0.56. 



(43) 



C. Dispersing point vison 

One can lower the energy by constructing plane waves. 
To compute the associated dispersion relation we have to 
evaluate the following matrix elements: 



m\j) 



We begin with the overlap matrix 



(44) 
(45) 



(46) 



(N(i) and N(j) are defined in Eq. 35) which simplifies to 



1 



% = ^]T< c |(-i) 



N(i,j) I 



TV ■ 

c 



(47) 
(48) 



where the local operator N(i,j) counts the number of 
dimers across some path fiy connecting the triangles i 
and j (see Fig. 2), and No(i,j) is equal to that num- 
ber in the reference configuration cq, chosen with all the 
dimers horizontal (Fig. 2). This follows from two simple 
properties: (c| ( — l) JV (' l )+ JV (j)+ JV (*.i) | c ) is independent of 

the configuration c), and (c |(-l)* ( ' i)+ ^ (j) |co) = 1. We 
finally write 



(-1) 



N{i,j] 



(49) 



where (• • •) represents the average with equal weight over 
all dimer coverings. 

The average of any such diagonal observable can be 
computed using the Pfaffian of a (modified) Kasteleyn 
matrix. 1,24 In the present case, the "string" observable 
(— is coded in the Kasteleyn matrix by changing 
the signs of the matrix elements corresponding to bonds 
crossed by f2.y, i.e., setting some bond fugacities to — 1. 
To evaluate numerically such an expectation value, we 
construct the modified Kasteleyn matrices corresponding 
to a large enough triangular lattice, in which the path 
f2ij is embedded. Finite-size effects decay exponentially 
with the system size, so that lattices with 28 x 28 sites 
(with periodic boundary conditions) can safely be used 
to evaluate Sfj with high accuracy up to distances d ~ 10 
between triangles i and j. 

The matrix elements S^ , plotted in Fig. 7 as a function 
of the distance between i and j, decay exponentially (a 
result anticipated by Read and Chakraborty~' J ). These 
quantities have already been evaluated by Ioselevich et 
al. 5 using a classical Monte Carlo sampling. 

By construction, a product like ^^^M ' ' ' Q„,ii 
(where • ■ ■ ,i n form a closed loop of triangles) is 

equal to the parity of the number of sites enclosed in 
the loop. Thus, the signs of the matrix elements 
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FIG. 7: (Color online) Overlap \ {i\j) \ between two point-vison 
states (denned in Eq. 39) as a function the distance dij be- 
tween i and j. The (blue) square correspond to Eq. 44 and the 
three (red) crosses correspond to Eq. 61 with a — —0.8. The 
calculations are done using an exact evaluation of the Pfaf- 
fians on a 28-site lattice with periodic boundary conditions. 
Dashed line: guide to the eye corresponding to an exponential 
decay with the dimer-dimer correlation length £ = 0.76. 
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arc similar to those of the hopping amplitude of a parti- 
cle moving on the hexagonal lattice and subjected to a 
magnetic field corresponding to half a flux quantum per 
hexagon. In such a case, the magnetic unit cell has to 
be doubled compared to the original lattice cell. Since 
the original unit cell contains one triangular site and two 
triangular plaquettes, the magnetic one contains four tri- 
angular plaquettes and thus four sites of the hexagonal 
lattice. In Fourier space, S°(k) is a 4 x 4 matrix, as the 
matrix M discussed previously. The same is also true for 
the Hamiltonian matrix elements described below. 

The spectrum of the overlap matrix S°(k) is plotted 
in Fig. 8 (for k describing a representative path in the 
Brillouin zone). As an important result, the eigenvalues 
are strictly positive for all k. To our knowledge, it is 
the first time that the linear independence of point vison 
states is proved. 

Let us now turn to the matrix elements of the Hamil- 
tonian. We wish to transform Eq. 45 into an expression 
which can be evaluated using the Pfaffians, that is the 
expectation value of a diagonal observable in the dimer 
basis. First, we write H. as a sum of projectors 



H 

K 



2 E^o 

ro 

IVVoXVVo 



V: 



(50) 

(51) 
(52) 



FIG. 8: (Color online) Spectrum of the overlap matrix S a (k) 
(Eqs.44 and 61), for different values of a (0, -0.2, -0.5 and 
-0.8), along the path A^B^C^A in the Brillouin zone (see 
Fig. 6). The bottom panel is a zoom on the lowest eigenvalues 
of S a . In these calculations, the matrix elements S"j are ne- 
glected for triangles i and j at distance d > d max = 10 (102 th 
neighbor on the hexagonal lattice). Up to this distance, the 
finite-size lattice used in the calculation (28 x 28 sites) gives 
practically the infinite volume limit for S"j . The value of d ma x 
used here is large enough to ensure a good convergence of the 
spectrum since the curves obtained with a smaller truncation 
distance (d max ~ 8.47 - shown here with the same colors) 
are almost superposed with that for d max = 10, except for 
the bottom of the spectrum at a = —0.5 and a = —0.8. 
The overlap spectrum turns out to be gapped for a = and 
a — —0.2 (and probably at a = —0.5 and a = —0.8 too), 
indicating the linear independence of the vison states. 



and expand Eq. 45 into 

2 EE(- 1 ) Ar(ci,l)+JV(C2j) ( c i|n I 'ok 2 ),(53) 



m 



ci,c 2 r 



with the notation N(c,i) = (c\N(i)\c) = ±1. Because of 
the double sum ^2 , this does not yet have the form 
of a diagonal observable amenable to an evaluation with 
Pfaffians. To go further we note that (c\ |II ro |ca) vanishes 
if ci and/or C2 is not flippable around the rhombus ro- 
In addition, (ci|II ro |c2) = if c\ and C2 differ anywhere 
outside r . So we can restrict the sum to pairs of con- 
figurations (c, c) which differ by a dimer flip in ro and 
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which are identical elsewhere on the lattice: 
H ij = JfYl E 



{<■■<■) 

F ro (c) = l 



n 



(-l) N ^ l \c\ + (-l) N ^(c\ 
(-l) N ^\c) + {-l) N ^ j 



(54) 



where F ro (c) = 1 if c is flippable on rhombus ro, and 
F ro (c) = otherwise. 

c and c only differ inside r , so the signs ( — l) N ( c ' 1 ) 
and (— 1) JV ( C ' 1 ) are the same if i is not inside ro, and are 
opposite if i e rp. Let us note 5,- iro = — 1 if i £ ro and 



^i,r — 1 otherwise. This leads to 



Hi = ^ ^ 



2 



(_l)^(c,i)+AT(cJ) 

(c,c) 

^r (c) = 1 

x[(c|+^,ro(5|]n ro [|c)+5 J>D |c)]. (55) 



Using the explicit form of li ro , we get 



[(c\+S t ^ ro (c\]U ro 



1 



c ) + 6 3,r„\c)} = ^(1 - <Vo)(l ~ 5 J>o) 



(56) 



and, finally, 
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FIG. 9: (Color online) Variational vison (Eq. 60) dispersion 
relation for different values of a, along the path A— +B— >C— - >A 
in the Brillouin zone (see Fig. 6). The curve for a = 
corresponds to the variational energies of point-vison states 
(Eq. 39). The absolute minimum is -E m i n = 0.119, found at 
(k x ,k y ) = (7r/6,7r/2) = B for a value ~ —0.8 of the varia- 
tional parameter a. The exact value of the gap (0.089) 1 ' is 
marked as a dotted horizontal line. In this calculation, the 
matrix elements S?j are neglected for triangles i and j at dis- 
tance d > dmax = 10. This value is large enough to ensure 
a perfect convergence of the spectrum below E < 0.6, as can 
be checked from the fact that the dispersion curves obtained 
with a smaller truncation distance (d max ~ 8.47 - shown here 
with the same colors) are practically identical. Some higher 
energy states, however, are not fully converged. 



' 2Af 



E 

Fr (c) = 1 



("I) 



N(c,i)+N(c,j) 



(57) 



1 



X 



2J\f 



(1 - S itro )(l - Sj. ra ) 

(c\(-lf^\c). (58) 



c 

F ro (c) 



The last expression is the average of a diagonal observable 
and can, thus, be evaluated using Pfaffians: 



(1 - ^i,ro)(l - Sj jTo ) 



X ((-!)*<■*•?> F r , 



(59) 



where we used an operator notation for the "flippabil- 
ity" F ro \c) — F ro (c)\c). To get a nonzero contribution, 
there must be at least one rhombus ro containing i and 
j. So H®j = if i and j are not nearest neighbors. Us- 
ing modified Kasteleyn matrices in a way similar to that 
described for the overlap matrix S°, the nonzero matrix 



elements can be calculated. We find \Hfj\ = 6po for i = j, 
\H®j\ = 2p when i and j are first neighbors. 

Although the matrix elements of TC are simple, the 
band structure is, however, not that of a simple tight- 
binding Hamiltonian, because of the non-orthogonality 
of the present variational vison states. In Fourier space, 
S°(k) and H°(k) are 4x4 matrices. The spectrum is 
obtained by solving the generalized eigenvalue problem 
H°(k)tp k - E(k)S°{k)ip k = 0. The results are shown 
in Fig. 9 (a = curve). The minimum of e k is found at 
(k x , k y ) = (tt/6, 7r/2), in agreement with the Monte Carlo 
results. 6 ' 18 The corresponding energy is e m ; n = 0.16, 
which is significantly lower than the energy (0.56) found 
above for a completely localized point vison \i). The bot- 
tom of this variational point vison band is, however, still 
high compared to Ivanov's 1 ' estimate (0.089) of the gap in 
the vison sector. The next section provides an improved 
family of variational states. 



D. Dressed vison 

The vison wave function of Eq. 39, differs from the UK 
wave functions only through minus signs. It has the "ac- 
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cidental" property that dimer-dimer correlations are the 
same as in the ground state." As a consequence, the exci- 
tation energy of such a state is carried only by the kinetic 
part of the Hamiltonian. Clearly, some local reweighting 
of the dimer configurations in the vicinity of the vortex 
core would allow an optimized balance between the po- 
tential and kinetic costs of the excitation, and would lead 
to an improved variational wave function. This amounts 
to "dressing" locally the initial point-vison state by some 
even - but fluctuating - number of additional point vi- 
sons. 

As a simple improvement of the vison wave function, 
we introduce a variational parameter a to reweight the 
configurations depending on their "flipability" at the core 
of the vison. This gives the following vison state: 

\i,a) = -^£(-l)*W(i + (60) 

where Fj | c) = | c) if the dimer configuration c is flippable 
around one of the three rhombi r±(i), r 2 (j), and r^{i) 
containing the triangle i, and i^|c) =0 otherwise. 

As for the point vison of Eq. 39, the vison states of 
Eq. 60 are not orthogonal and we have to evaluate their 
overlaps: 

Sg = <i,a|i,a) (61) 
Repeating the transformation leading to Eq. 49, we get 

^ = ^•((-1)^(1 + 0^) (l + aFj)), (62) 

which is an expectation value for a diagonal operator 
that we evaluate using Pfaffians. In addition to the sign 
changes due to (— 1) JV ( , 'J), some entries of the Kasteleyn 
matrix have to be modified to incorporate the flippability 
operators F r . More precisely, counting only the coverings 
which satisfy F r = 1 is done by "isolating" the rhombus 
r, that is, by switching to zero in the Kasteleyn matrix 
the 14 bonds which connect the sites of rhombus r to 
their neighbors outside r. 

Beyond some distance between the vison cores i and j 
(fourth neighbor on the hexagonal lattice), no rhombus 
can touch simultaneously both triangles. In that case, it 
can be shown that Sfj — Sfj is independent of a. The 
eigenvalues of the overlap matrix S a (k) are displayed in 
Fig. 8 for a few selected values of a. Although a full con- 
vergence as a function of the truncation distance d max 
(see caption of Fig. 8) has not been obtained, we believe 
that there is a finite gap for all the values of a shown 
here, and that the dressed vison states are linearly inde- 
pendent. 

The evaluation of the Hamiltonian matrix elements 

H^ = (i,a\H\j,a) (63) 

for dressed vison can still be done using Pfaffians, but the 
algebraic manipulations are slightly more lengthy than 



the previous ones and we refer the reader to Appendix B. 
We find nonzero matrix elements H°j up to (and includ- 
ing) the 9 th neighbor (compared to first neighbor for 
point vison). 

The results for the variational dispersion relation are 
shown Fig. 9. The qualitative shape of the lowest band 
is almost unchanged compared to the point vison states 
(a = 0), except for an almost uniform shift which lowers 
the gap. The minimum of Ek is found at (k x ,k y ) = 
(7r/6,7r/2) = B for a variational parameter a ~ —0.8, 
and the corresponding energy is E min = 0.119, about 
33% higher than the exact value. 

The fact that simple wave functions like those of Eq. 39 
or Eq. 60 reproduce qualitatively the shape of the exact 
dispersion relation is presumably due to the fact that the 
dimer-dimer correlation length is rather small at the RK 
point (of order of one lattice spacing 1 ). As a consequence, 
the exact vison states only differ from Eq. 39 at short dis- 
tances from the vortex core, and the long-distance part 
(string), responsible for the flux tt per hexagon, is essen- 
tially exact. This is, of course, no longer true away from 
the RK point and in the direction of the crystal, where 
the correlation length rapidly grows. 

In an attempt to extend this variational approxima- 
tion away from the RK point, we computed the gap of 
the dressed vison states in perturbation theory to first 
order in (1 — V). However, at this order, the gap turns 
out to close very slowly and does not lead to a mean- 
ingful estimate for the critical V at the liquid/crystal 
transition. This failure is closely related to the remark 
above: the size of the core of the vison presumably grows 
rapidly away from the RK point, a feature which cannot 
be accounted for with the present variational states. 



IV. CONCLUSIONS 

In this paper, we have developed two simple ap- 
proaches to describe the vison excitations of the QDM 
on the triangular lattice. The first one is based on 
a soft-dimer version of the model, which exactly takes 
the form of a 1^ gauge theory. We have shown that a 
semiclassical spin-wave approximation to the spectrum 
of the dual Ising theory captures the important fact that 
the disappearance of the RVB liquid is due to a vison 
condensation. It also reproduces the qualitative shape 
of the vison spectrum in the disordered phase and, more 
importantly, at the transition (linear spectrum at the cor- 
rect points of the Brillouin zone), as well as the spatial 
pattern of the ordered crystalline state. 

The second approach is a variational approximation to 
the vison wave functions at the RK point of the (hard- 
core) QDM. It reproduces semiquantitatively the vison 
dispersion relation, and provides a simple picture for the 
vison wave functions which goes beyond the naive point- 
vison approximation. 

Beyond the problem of vison dispersion and condensa- 
tion, we expect these approaches to be useful for other 
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problems related to vison excitations in QDM, in partic- 
ular, their mutual interaction and their interaction with 
vacancies and mobile holes. This is left for future inves- 
tigations. 
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APPENDIX A: FOURTH ORDER EFFECTIVE 
QUANTUM DIMER MODEL 

The goal of this Appendix is to derive an effective 
Hamiltonian for the model of section II. A defined by the 
Hamiltonian 



H = Hj + H r = -jJ2rf-Tj2U T ii 



(Al) 



Hi) 



in the limit Y / J -C 1 to fourth order in T/ J. For simplic- 
ity, let us define the unperturbed Hamiltonian Hq = Hj 
and the perturbation V = Hp. Since the Hilbert space of 
the model is restricted by the constraint that the number 
of dimers starting from a given site is odd, the ground 
state manifold of the unperturbed Hamiltonian for pos- 
itive J consists of all states having exactly one dimer 
emanating from each site. Let us denote by Pq the pro- 
jector onto the Hilbert space generated by these states, 
and by S the resolvent defined by 



S=- 



Ho — Eq 



(A2) 



It is easy to check that V changes the parity of the num- 
ber of dimers. Indeed, the term of V acting on the tri- 
angle i transforms the states with and 3 (resp. 1 and 
2) dimers around the triangle i into each other. This 
implies that PqVPo = 0, and more generally that the 
effective Hamiltonian only contains even powers of the 
perturbation. Thanks to the property PqVPq = 0, the 
fourth order contribution reduces to 3 terms, and the 
effective Hamiltonian up to fourth order reads: 

H cS = PqVSVPq + PqVSVSVSVPq 

~ 2 P o (VS 2 VP VSV + VSVP VS 2 V) P 

Up to a constant, this effective Hamiltonian can be 
written as a QDM acting on 4 and 6-site plaquettes. The 
4-site Hamiltonian has the form of the regular RK model 






FIG. 10: The three types of 6-site plaquettes around which 
dimer shifts are generated at fourth order in T/J. 

with t = T 2 /J - T 4 /J 3 and V = T 4 /2J 3 . The 6-site 
Hamiltonian only consists of kinetic terms that flip the 
dimers around the three possible types of 6-site plaque- 
ttes shown in Fig. 10, with amplitude — 3r 4 /4J 3 for types 
(1) and (2), and with amplitude — r 4 /J 3 for type (3). 



APPENDIX B: HOPPING AMPLITUDE FOR 
THE DRESSED VISONS 

For the dressed vison states (Eq. 60) , the equivalent of 
Eq. 55 is 



'■" (c,c) 
F ro (c) = l 



N(c,i)+N(c,j) 



n,. 



c|(l + a/ , i)+<5i,r <c|(l + a£i) 
(1 + aF 3 )\c) + S 3 ,, (l + aF 3 )\c) (Bl) 



Using the explicit form of the projector 1%, we get 
1 

= 2A7 



(i,a\IL ro \j,a) 



(_l)Ar(c,i)+JV(c,i) 

(c, c) 

Fr (c) = l 

x (1 + aFi{c) - 5 i>ro {l + aFi{c))) 
x(l + a F J -( C )-^ rD (l + aF i (c))|B2) 



Unlike Eq. 57, both the coverings c and c enter the 
expression. This does not have the form of a di- 
agonal observable and the terms F r (c) (with r € 
{ri(i),r 2 (i),r 3 (i),r 1 (j),r2(j),r 3 (j)}) need to be elimi- 
nated to allow for an evaluation with Pfaffians. By in- 
specting the possible relative positions of two rhombi r$ 
and r, one arrives at the following two relations: 

• If two rhombi ro and r have 0,1,3 or 4 sites in com- 
mon, F r (c) = F r (c) for any pair of configurations 
(c, c) which differ by a flip around the rhombus ro- 

• It ro and r have 2 sites in common, let us call b the 
bond of r which does not touch ro. Then we have 
F r (c) = D b (c)(l - F r (c)), where D b (c) = 1 if b is 
occupied by a dimer of c, and otherwise. 
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We may combine the two cases above into some com- 
pact notation F r (c) = A ryTQ {c), valid for any pair of 
configurations (c, c) which differ by a flip around the 
rhombus r . Accordingly, we define an operator A i ro — 

^7-1(1)^0+^2(1)^0+^3(1)^0 f° r each triangle i and rhom- 
bus r . The matrix element of Eq. Bl is now expressed 
as the expectation value of a diagonal operator: 

h?j = |E((- 1 )* (,J) ^ 

x (l + aFi - S i<ro (l + aAi iro )\ 

X (l + aFj - 5 jtro (l + aA j>ro )\ ) (B3) 

This expression has to be expanded into a polynomial 
in F and D operators before each term can be evalu- 



ated thanks to the Pfaffian of an appropriate Kasteleyn 
matrix. As before the (— l) Ar ( 1 ^) introduces some sign 
changes, and each F (or D) operator requires to "iso- 
lating" the corresponding rhombus (or bond) by switch- 
ing to zero the corresponding matrix elements. Several 
tens of terms typically appear for each pair of triangles 
{if), and an automated treatment by computer had to be 
coded to obtain the hopping amplitudes. The results of 
Fig. 9 represent several hundreds of CPU hours using the 
software Maple to generate all the correlators and eval- 
uate the corresponding Pfaffians on a finite lattice with 
28 x 28 sites. 
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|(1 — a x ) (as here for visons), magnetic long-range order 
in the z (or y) direction (for the spin variables) is equiva- 
lent to Bose condensation (off-diagonal long-range order in 
the Boson operators). An important difference is, however, 
that the number of visons is not conserved. Only their par- 
ity is conserved. Accordingly, the vison condensed phase 
spontaneously breaks a discrete (Z2) gauge symmetry (<r z , 
which is not gauge invariant, acquires a finite expectation 
value in the crystal phase), and not a continuous [U(l)] 
one. Consequently, the condensed phase is gapped and does 
not have a Goldstone mode. In the dimer model, the spec- 
trum is, therefore, gapless only at the transition. 



